#1、包的安装
#输入命令：
>source(“http://bioconductor.org/bioLite.R”)
>bioLite(“DESeq”)
#下载安装相应的软件包。

#2、数据
>library("DESeq")    		#加载DESeq包
>Sample_set=read.table(file ="sample_set.txt",header=TRUE)			#读入准备好的数据
> condition=colnames(Sample_set[-1])
> condition			#得到样本名
#[1] "GFP"      "PsCRN115" "PsCRN161"
> cds=newCountDataSet(Sample_set[-1],condition)			#将数据利用函数读入变量cds，新的数据结构
#数据标准化
> cds=estimateSizeFactors(cds)			#函数estimateSizeFactors估计统计数据的大小因子
> sizeFactors(cds)
#      GFP  PsCRN115  PsCRN161 
#0.9079696 1.0000000 1.1226796 
> head(counts(cds,normalized=TRUE))			#如果统计数据的每列除以这列的大小因子，这样统计值就变成同一规模，使它们具有可比性。函数counts可以做这个计算。
#          GFP PsCRN115   PsCRN161
#1   13.216302        6   8.907261
#2   46.257058       20  26.721783
#3    1.101359        3   4.453630
#4 1083.736778     1067 660.028031
#5    9.912227       13   9.797987
#6  148.683399      231 144.297626

#3、方差估计
#DESeq推断依靠估计典型的数据间的方差和平均关系，或者等效的离差和均值。离差可以理解为生物变异系数的平方。估计离差可以使用以下命令。
> cds=estimateDispersions(cds,method="blind",sharingMode="fit-only")			#函数estimateDispersions做了三步，首先估计每条基因的离差，然后，通过估计匹配一条曲线，最后，每个基因分配一个离差，从每条基因估计值和匹配值选一个。为了让用户知道中间过程，fitInfo对象被储存下来。
>str(fitInfo(cds))
#画图并保存，函数plotDispEsts可以画出每条基因的估计值和平均正常统计值的关系。
> pdf('cds.pdf')
> plotDispEsts(cds)
> dev.off()

> head(fData(cds))
  disp_blind
1 0.19444528
2 0.16149828
3 0.49344090
4 0.06503826
5 0.17611634
6 0.07248443

#4、差异表达基因，两个样本的标准比较
#GFP/PsCRN115
>GP115=nbinomTest(cds,"GFP","PsCRN115")
>GP115=data.frame(Sample_set$gene_names,GP115)			#添加gene_names
>head(GP115)
#我们先画log2折叠变换和平均正常统计量的关系，红点表示在10%FDR的基因。
>pdf('GP115.pdf')
>plotMA(GP115)
>dev.off()
 
#根据P值统计直方图
>pdf('GP115_hist.pdf')
>hist(GP115$pval,breaks=100,col="skyblue",border="slateblue",main="")
>dev.off()
 
#通过FDR过滤掉基因
>GP115_Sig=GP115[GP115$padj<0.1,]
#列举最有效的差异表达基因
>GP115_differential=GP115_Sig[order(GP115_Sig$pval),]
>head(GP115_differential)
#也可以看到有效基因的最强下调
>GP115_down=GP115_Sig[order(GP115_Sig$foldChange,-GP115_Sig$baseMean),]
>head(GP115_down)
#最强上调
>GP115_up=GP115_Sig[order(-GP115_Sig$foldChange,-GP115_Sig$baseMean),]
>head(GP115_up)
#保存输出到指定CSV文件
write.csv(GP115_differential,file="GP115_differential.csv")
#并标记上调下调基因，得到GP115_differential.xlsx。

#GFP/PsCRN161
>GP161=nbinomTest(cds,"GFP","PsCRN161")
>GP161=data.frame(Sample_set$gene_names,GP161)
>head(GP161)
#我们先画log2折叠变换和平均正常统计量的关系，红点表示在10%FDR的基因。
>pdf('GP161.pdf')
>plotMA(GP161)
>dev.off()

#根据P值统计直方图
>pdf('GP161_hist.pdf')
>hist(GP161$pval,breaks=100,col="skyblue",border="slateblue",main="")
>dev.off()

#通过FDR过滤掉基因
>GP161_Sig=GP161[GP161$padj<0.1,]
#列举最有效的差异表达基因
>GP161_differential=GP161_Sig[order(GP161_Sig$pval),]
>head(GP161_differential)
#也可以看到有效基因的最强下调
>GP161_down=GP161_Sig[order(GP161_Sig$foldChange,-GP161_Sig$baseMean),]
>head(GP161_down)
#最强上调
>GP161_up=GP161_Sig[order(-GP161_Sig$foldChange,-GP161_Sig$baseMean),]
>head(GP161_up)
#保存输出到指定CSV文件
write.csv(GP161_differential,file="GP161_differential.csv")
#发现样本GFP与样本PsCRN161之间没有差异表达基因。

#单一样本不用这一步：
#对于弱表达基因，需要加更强的改变，这边以两个样本基因为例。
> ncu1= counts( cds, normalized=TRUE )[ , -3]			#ncu1为GFP、PsCRN115两列矩阵
>head(ncu1)
>pdf(‘GP115_mean.pdf’)
>plotMA(data.frame(baseMean = rowMeans(ncu1), + log2FoldChange = log2( ncu1[,2] / ncu1[,1] )), +  col = "black")
>dev.off()

> ncu2= counts( cds, normalized=TRUE )[ , -2]			#ncu2为GFP、PsCRN161两列矩阵
>head(ncu2)
>pdf(‘GP161_mean.pdf’)
>plotMA(data.frame(baseMean = rowMeans(ncu2), + log2FoldChange = log2( ncu2[,2] / ncu2[,1] )), +  col = "black")
>dev.off()